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Recently in functional magnetic resonance imaging (fMRI) studies there has been 
an increased interest in understanding the dynamic manner in which brain regions 
communicate with one another, as subjects perform a set of experimental tasks or as 
their psychological state changes. Dynamic Connectivity Regression (DCR) is a data-driven 
technique used for detecting temporal change points in functional connectivity between 
brain regions where the number and location of the change points are unknown a priori. 
After finding the change points, DCR estimates a graph or set of relationships between 
the brain regions for data that falls between pairs of change points. In previous work, 
the method was predominantly validated using multi-subject data. In this paper, we 
concentrate on single-subject data and introduce a new DCR algorithm. The new algorithm 
increases accuracy for individual subject data with a small number of observations and 
reduces the number of false positives in the estimated undirected graphs. We also 
introduce a new Likelihood Ratio test for comparing sparse graphs across (or within) 
subjects; thus allowing us to determine whether data should be combined across 
subjects. We perform an extensive simulation analysis on vector autoregression (VAR) 
data as well as to an fMRI data set from a study (n = 23) of a state anxiety induction 
using a socially evaluative threat challenge. The focus on single-subject data allows us to 
study the variation between individuals and may provide us with a deeper knowledge of 
the workings of the brain. 
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1. INTRODUCTION 

Traditionally, functional magnetic resonance imaging (fMRI) 
studies have concentrated on locating brain regions showing 
task-related changes in neural activity, for example, greater activ- 
ity during an experimental task than during a baseline state 
(Lindquist, 2008). The voxel-wise general linear model (GLM) 
(Worsley and Friston, 1998) has become the standard approach 
for analyzing such data. However, in order to gain a deeper insight 
into the neural processing of the brain as the subject performs a 
set of experimental tasks or changes psychological state, there is 
now an intense interest in the field to study the interactions or 
connectivity between distinct brain regions. 

A very active area of neuroimaging research involves exam- 
ining the directed (effective connectivity) and undirected asso- 
ciation (functional connectivity), between two or more spatially 
remote brain regions (Friston et al., 1993; Biswal et al, 1995). For 
example, the "Human Connectome" project is currently uncov- 
ering the structural and functional determinants of brain activity 
using information from anatomical, functional and effective con- 
nectivity (Sporns et al., 2005). For fMRI data, the covariance or 
correlation matrix as well as the precision matrix (inverse covari- 
ance matrix) can be estimated and used to measure functional 
connectivity between different regions based on their time series. 



However, these techniques are only useful when estimating static 
functional connectivity networks or when the experimental con- 
dition does not change over time. The techniques do not consider 
the possibility of changing functional connectivity or the dynam- 
ical nature of interactions between brain regions. If the dynamic 
connectivity (or dependencies) between brain regions are not 
taken into account, the true dynamics are not captured and the 
overall connectivity structure of the data will simply represent 
an aggregate of all the changing connectivities between runs. 
Specifically, the distinction between the different psychological 
states would not be found, rather an average of the functional 
connectivity would be estimated, which in many experiments 
would lead to meaningless results. 

The popular PPI (Psychophysiological Interactions; Friston 
et al., 1997) method addresses the dynamic nature of functional 
connectivity but it assumes the timing of the various contexts or 
state changes is known. However, often, it is difficult to specif)' 
the nature, timing and duration of the psychological processes 
being studied a priori. Chang and Glover (2010) and Allen et al. 
(2013) introduced methods for estimating time varying connec- 
tivity between brain regions for resting state data using a moving 
window. Hence, an important extension of these works is to 
have methods that can detect changes in connectivity, regardless 
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of the nature of the experimental design. In order to find evi- 
dence of fluctuations in functional connectivity over the course 
of the experiment, we previously introduced a technique called 
Dynamic Connectivity Regression (DCR; Cribben, 2012; Cribben 
et al., 2012). DCR is a data-driven method for detecting func- 
tional connectivity change points between brain regions where 
the number of change points and their location are not known 
in advance. After finding the change points, the method estimates 
a graph or series of relationships between the brain regions for 
data in the temporal partition that falls between pairs of adja- 
cent change points. It is assumed that the graph does not change 
within each partition. The motivation behind the method stems 
from the fact that under certain circumstances the dynamic man- 
ner in which the brain regions interact with one another over time 
is of primary interest. DCR utilizes the graphical lasso or glasso 
(Friedman et al., 2007) and so it allows for change point and 
graph estimation in a high-dimensional context by setting many 
elements of the precision matrix to zero. DCR can be applied to 
estimate the dynamic interaction between the brain regions as 
well as the timing, duration and details of the interaction. 

In Cribben et al. (2012), DCR was used to estimate change 
points and undirected graphs for multi-subject data. In this paper, 
we focus exclusively on single-subject data because group aver- 
age analyses risk removing unique individual patterns of activity 
and analyzing single-subject data allows us to understand the 
sources of intersubject variability in brain activity. It is assumed 
that by averaging across the group, the noise effects are removed 
but by analyzing individual subjects and the variations between 
them, the results may provide insight to the understanding of 
the mind/brain relationship. In other words, while analyses of 
multi-subject data are highly valuable for understanding general 
cognitive processes, the study of variation between individu- 
als may be able to provide us with a deeper knowledge of the 
workings of the brain. Furthermore, by analyzing multi-subject 
data we implicitly average connectivity patterns across the group, 
which may lead to meaningless results. An example of this can 
be found in section 5, where subjects with positive connectivity 
is canceled out by subjects with negative connectivity. Evaluating 
how functional brain networks change over time and relating this 
information to brain dysfunction promises to be very important 
for understanding the underlying mechanisms of brain disorders 
such as depression and Alzheimer's disease. 

To this end, we introduce a new DCR algorithm that involves 
the addition of extra steps to the old algorithm to increase accu- 
racy for single-subject data with a small number of time points. 
The extra steps prove vital in situations where there are small and 
subtle connectivity changes. We show that the new DCR algo- 
rithm is very adept at finding both major and minor changes 
in the functional connectivity structure. A comparative study 
between the old and new algorithm is carried out in the simu- 
lation study. The stationary bootstrap inferential procedure for 
detecting significant change points is also revised accordingly to 
obtain more accurate estimates of the distribution of the BIC 
reduction at each candidate change point. The glasso technique 
for estimating sparse graphs can be severely impacted by the 
inclusion of only a few contaminated values, such as spikes that 
commonly occur in fMRI time series, and the resulting graph 



has the potential to contain false positive edges (Finegold and 
Drton, 2011). The second extension involves performing a boot- 
strap inferential procedure on the edges or partial correlations in 
the undirected graphs to remedy the issue of the false positive 
edges. 

By carrying out the analysis on individual subjects we do not 
need to collapse information across subjects or assume that the 
graph structure between each pair of change points is the same 
across subjects. The subject's change points and partition specific 
undirected graphs are thereafter compared to determine whether 
the change points and connectivity structures are common across 
all subjects. We also introduce a new Likelihood Ratio test on 
the similarity of the sparse precision matrices across subjects. 
With this information we can make inferences about whether 
the assumption of equal connectivity structure across subjects 
is valid. The matter of whether to combine information across 
subjects (to analyse group average results) or whether to look 
at individual patterns of brain activity is critical for individual 
treatment of brain disorders and very topical in the field of neu- 
roscience (see Van Horn et al, 2008 for a discussion). As noted in 
Varoquaux et al. (2010), the use of multiple datasets is very chal- 
lenging due to the variability between subjects coupled with the 
variability of functional signals between experimental runs. 

This paper is organized as follows. In section 2 we begin by 
briefly discussing the theoretical foundations of DCR. Next we 
describe the new DCR algorithm for single-subject data with a 
small number of time points and inferential procedures for deter- 
mining statistically significant connectivity change points and 
statistically significant edges in the undirected graphs. We also 
discuss a Likelihood Ratio test for comparing precision matri- 
ces across and within subjects. In sections 3 and 4 we outline the 
simulations and the fMRI study that the new DCR is applied to 
respectively. We describe the results of the simulations and exper- 
imental data in section 5. The results display DCR's capability to 
find connectivity change points for single-subject data with only 
a limited number of observations thus allowing for comparison 
across subjects. 

2. MATERIALS AND METHODS 

2.1. DYNAMIC CONNECTIVITY REGRESSION (DCR) 

A more detailed explanation of the theory behind DCR can be 
found in Cribben et al. (2012) and in Cribben (2012). Suppose 
Y (e.g., several brain regions) has n i.i.d data points sampled 
from a zero mean p-dimensional multivariate Gaussian distri- 
bution, iV(0, E). In order to estimate the dependence structure 
of Y we could use the sample covariance matrix, E = (« — 
1) ' X!i'=iO',- ~ j')0'/ ~ y)^ (estimate of E). Alternatively we 
could use the precision matrix (inverse covariance matrix), Q, = 
E^^ The natural estimator for Q is S^'. However, we could 
also consider regularization and sparse techniques which are espe- 
cially needed when p > n. The most widely used precision matrix 
estimation procedures are based on penalized likelihood maxi- 
mization or constrained error minimization approaches which 
are optimization algorithms with different objective functions. 
All of these estimation procedures assume that the dependence 
structure between the data remains constant. However, in many 
applications, including neuroimaging, the dependency structure 
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between the individual voxel or ROI time series varies over time. 
If the dynamic dependencies are not properly accounted for, 
the overall dependence structure will simply be an aggregate 
of all the different dependency structures. DCR detects change 
points in functional connectivity between brain regions, without 
knowledge on the number or location of the change points, and 
estimates a graph or set of relationships between the brain regions 
for data between each pair of change points. 

DCR builds upon the extensive literature on graphical mod- 
els (Whittaker, 1990; Edwards, 1995; Cox and Wermuth, 1996). 
Within this framework, graphical models display the depen- 
dency structure of a random vector Y using a graph G. Graphs 
are mathematical structures that can be used to model pair- 
wise relationships between variables. They consist of a set of 
vertices V and corresponding edges E connecting pairs of ver- 
tices. Here, we focus exclusively on undirected graphs which do 
not infer the directionality of dependence or functional con- 
nectivity between the variables or brain regions. A graph of 
Y can alternatively be represented using the precision matrix 
(inverse covariance matrix) of Y, with the elements of the matrix 
corresponding to edge weights. Here a missing edge between 
two vertices in the graph indicates conditional independence 
between brain regions, giving rise to a zero element in the preci- 
sion matrix. Throughout this paper we will model dependencies 
between brain regions using the precision matrix or undirected 
graphs and we wOl use them interchangeably throughout the 
paper. 

2.2. THE NEW DCR ALGORITHM 

In this work, we introduce a new DCR algorithm necessary for the 
efficient and accurate analysis of single-subject fMRI data. The 
step-by-step description of the new algorithm can be found in 
section A.l. of the Appendix. The new steps involve recalculat- 
ing the Bayesian Information Criterion (BIG: Equation 4) at each 
candidate change point after the first pass of the algorithm, alter- 
ing the stationary bootstrap inference procedure on the candidate 
change points and including a new bootstrap inference proce- 
dure on the undirected graphs within each partition in order to 
remove the false postive edges prevalent in undirected graphs for 
single-subject data. 

The objective of the algorithm is to split the time course or the 
time domain, X, into partitions Xi>- ■ ■ , X.m- As before, our greedy 
procedure begins by calculating the sample mean and sample 
covariance matrix for the entire data set using Equations (1, 2), 
respectively. 



E,= i-f(^. e Xj) , = 1 



(1) 



(2) 



For the entire data set calculation x.; is all the time points and tj is 
the number of time points in the entire experiment. Estimates of 
the precision matrix, f2, for the entire data set are then obtained 



using the graphical lasso or glasso (Friedman et al., 2007) for 
a specified full path of X values using Equation (3). The li- 
penalty in this equation induces sparsity and regularization on 
the elements of the estimated precision matrix. 



n : 



argmin^^Q {tr - log|^2| -|- 



(3) 



The tuning parameter A controls the sparsity of the precision 
matrix, large values giving rise to a very sparse precision matrix 
and small values giving rise to a very "full" precision matrix or 
graph. An efficient algorithm has been developed (Friedman 
et al., 2007) for finding Q, an estimate of the precision matrix, 
that estimates a single row (and column) of Q in each iteration by 
solving a lasso regression. We chose the value of X that minimizes 
the BIG 



sic=I] [(tr {y, - i^) {y, + k-\og (f,) 

i = 1 

(4) 

where tj is the number of data points in the partition and k 
is the number of off-diagonal elements in the precision matrix. 
After estimating the precision matrix based on BIG minimization 
and identifying non-zero edges, the model is refit without the h- 
constraint while keeping the zero elements in the matrix fixed to 
reduce bias and improve the model selection performance (Liu 
etal, 2010). 

The minimum BIG score for the entire data set is recorded, 
providing a baseline in which to evaluate subsequent splits. Note 
that the estimated precision matrix for the entire data set is simply 
an aggregate of all the changing connnections (or partial correla- 
tions) assuming there exists a dynamic connectivity structure in 
the data. The BIG is a model selection criterion based on combin- 
ing the likelihood function with a penalty term that guards against 
over-fitting. Hence, it balances the dual needs of adequate model 
fit and model parsimony. 

Upon completion of this step, the data is then partitioned into 
two sets; a left subset consisting of time points {1:A} and a right 
data set consisting of { A -|- l:T}, where T represents the length of 
the time course. The choice of A is of particular importance as 
it also represents the minimum possible distance between adja- 
cent change points and its value can be adjusted depending on 
the existence of a priori knowledge about the spacing of changes 
in functional connectivity. Unlike the multi-subject analysis pre- 
sented in Cribben et al. (2012), where DCR could find change 
points with close proximity (i.e., with A = 10), the single-subject 
analysis in this paper requires A to be a minimum of 30-40 
time points to ensure that there is enough data to provide stable 
estimates of the sample means and sample covariance matrices 
and hence the precision matrices. The sample mean and covari- 
ance matrices for both subsets are calculated separately using 
Equations (1, 2), respectively. Again the full regularization path 
of k values for the glasso is run and the values of k, and corre- 
sponding precision matrices, are chosen based on minimizing the 
BIG. Each model is then refit, as above, keeping the zero elements 
in each precision matrix fixed, and the combined BIG scores for 
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the two subsets are recorded. This procedure is repeated along the 
entire time path, with the data partitioned into two subsets with 
spht points ranging from {1:A+1} to {1:7 — A}. The partition 
with the smallest combined BIC score is chosen and, if its value 
is less than the BIC score for the entire data set, the correspond- 
ing time point is used to partition the data and the first candidate 
change point is identified. 

The DCR procedure continues by recursively applying the 
same method to each individual partition element until they can 
no longer be split any further. In other words, if the first candidate 
change point occurs at time p, the procedure is repeated for both 
the data set consisting of time points {l:p} and the one consist- 
ing of time points {p -|- 1:T). The procedure is repeated until no 
further splits reduce the BIC score. 

Once the candidate change points have been found, a new 
extra refitting step is carried out to improve the estimation of 
the BIC reduction at each split. In other words, the first candi- 
date change point separates the entire data set into two similar 
data sets based on their functional connectivity structure. But 
if further candidate change points are found then the first can- 
didate change point's BIC reduction is inaccurate and requires 
reestimation. This is particularly important for single-subject 
experimental data sets with few time points. When the DCR algo- 
rithm recursive procedure is finished, the candidate change points 
are ordered in time. If the first candidate change point occurs at 
time point p and the second candidate change point occurs at 
time point v, the BIC is calculated for data sets {l:p}, {p -|- l:v} 
and {l:v} by estimating the sample means, sample covariance 
matrices and precision matrices along the full path of A values 
for each partition. The new BIC reduction at time point p is then 

BIC{1 : v} - [BIC{1 : p} -|- BIC{p -I- 1 : v}] (5) 

By sequencing through all the candidate change points in this 
manner and recalculating the BIC reduction the model's accuracy 
is improved (compare Figures 4A,B as well as Figures 5A,B). This 
extra set of calculations only marginally increases the algorithm's 
computation time. After completion of the algorithm, the new 
DCR will have split the experimental time course into connected 
partitions xi, • ■ • , Xm within each partition X;i the glasso is 
again used to estimate an undirected graph Gj. 

The idea behind the DCR algorithm involves comparing the 
BIC of the entire data set with the BIC of two subsets (right 
and left) of the data. In the former case the model contains k 
parameters while in the latter case the model contains fci -|- k2 
parameters, where ki and ki represent the number of param- 
eters in the left and right subset, respectively. This allows for 
more free parameters in the two subsets case. However, this 
bias is balanced out by the fact that the full data model con- 
tains T observations while the split data model contains some 
combination of data that sums to T. It is evident that the BIC 
penalizes the split data model more than the BIC of the entire 
data set as the goal is to minimize the BIC. Hence, false pos- 
itive BIC reductions or false positive change points are not of 
great concern to the DCR algorithm. For more details see Cribben 
(2012). 



2.3. DETECTING SIGNIFICANT CHANGE POINTS 

To determine whether significant change points exist, confi- 
dence bounds for each non-zero BIC reduction at each splitting 
time are created using a new stationary bootstrap procedure 
(Politis and Romano, 1994). The stationary bootstrap, unlike 
a permutation procedure, is mindful of the dependency struc- 
ture inherent in the data where successive time points are 
assumed to be correlated, but observations "far apart" uncor- 
related. The stationary bootstrap is an adaptation of the block 
bootstrap (Liu and Singh, 1992; resamples blocks of observa- 
tions instead of individual observations) that allows for randomly 
varying block sizes. The new stationary bootstrap procedure 
is similar to the stationary bootstrap procedure introduced in 
the original DCR paper (Cribben et al, 2012), however, it is 
altered to reflect the changes in the DCR algorithm described 
above. 

The new stationary bootstrap procedure begins by running 
the new DCR algorithm to determine aU the ordered candidate 
change points. Assume the first and second ordered candidate 
change points are p and v, respectively. To test whether the BIC 
reduction at change point p is significant, data within the parti- 
tion { 1 :v} is repeatedly resampled across time using the stationary 
bootstrap. For each resampling, the data is split into two parts, 
one consisting of time points {l:p} and the other of {p -|- l:v}. The 
combined BIC score from each subset is subtracted from the BIC 
of the entire resampled data set {l:v} using the stationary boot- 
strap and the results are combined across replications to create 
a stationary bootstrap distribution for the BIC reduction at the 
candidate change point p. The procedure is repeated for each non- 
zero BIC reduction or candidate change point. For example, in 
order to check significance for the second candidate change point 
V, data within the partiton {p -|- l:t}, where 6 is the third candi- 
date change point, is repeatedly resampled across time using the 
stationary bootstrap. 

The (1 — a/2) and a/2 quantUes of the stationary bootstrap 
distribution of each non-zero BIC reduction for each splitting 
time can be plotted and interpreted as 100(1 — a)% confidence 
bounds. For a candidate change point at time p, if the BIC 
reduction for the original data is more extreme (either larger 
or smaller) than the 100(1 — a)% confidence bounds computed 
using repeated stationary bootstrap replicates, we conclude there 
is a significant change point at p, indicating a change in con- 
nectivity. This new stationary bootstrap procedure decreases the 
computation time of the algorithm as the size of the resampled 
data sets are now smaller. 

2.4. CONSTRUCTING UNDIRECTED GRAPHS 

Once the significant change points have been found, the data is 
divided into partitions defined by the splits. For each partition, 
the sample mean and sample covariance matrbc are estimate using 
Equations (1, 2), respectively. Using the full path of X values in 
Equation (3), the optimal X value and precision matrix is chosen 
based on minimizing the BIC. It is necessary to carry out these 
extra calculations as the sample means, sample covariance matri- 
ces, precision matrices, BIC reductions and X values associated 
with each partition calculated using the new refitting scheme in 
the DCR algorithm assumes that all candidate change points are 
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significant. If not, the estimated values associated with the par- 
tition are no longer valid as the time points in certain partitions 
change as a result. 

In order to make the off-diagonal elements of the precision 
matrix, P = Q, more interpretable they are converted into partial 
correlations, r,-,, using the following formula 



The undirected graphs in this work plot the partial correlations 
between the brain regions. 

2.5. TESTING EDGES 

A shortcoming of the glasso is the excessive number of false posi- 
tive edges present in the estimated undirected graphs. This fact is 
highlighted for multi-subject data in Cribben et al. (2012), but it 
is noteworthy that the false positive partial correlations tend to be 
small. The number of false positive edges increases in the single- 
subject setting, which is the emphasis here, with an excess of 
edges also complicating interpretation of the undirected graphs. 
Hence, we include an additional step in the new DCR algorithm 
by performing inference on the edges using information from a 
bootstrap procedure. 

The procedure is carried out for each partition specific undi- 
rected graph, computed as described in section 2.4. For a partition 
Xj, a bootstrap data matrix Y* . is formed by randomly selecting 
with replacement rows of data from the original partition specific 
data matrix Y^^.. We estimate a precision matrix (or undirected 
graph) for the new data set Y*^. for the full path of X values and 
we choose X and the undirected graph based on minimizing the 
BIC. We record each partial correlation in the undirected graph. 
A graph containing 20 vertices has 190 (0.5 x p x (p — I)) pos- 
sible edges or partial correlations. This procedure is repeated a 
large number of times (say, 1000). We then count the number of 
times glasso estimates the edge to be non-zero in the 1000 boot- 
strap runs. These proportions represent the probability of edge 
selection. For example, if glasso estimates the edge to be non-zero 
in 990 bootstrap runs, the edge has an estimated selection prob- 
ability of 0.99. For each edge, if its selection probability is greater 
than some pre-defined threshold, say 0.95, then the edge or partial 
correlation is significant and it remains in the undirected graph. 
On the other hand, if its estimated selection probability is less 
than the pre-defined threshold then the edge or partial correlation 
is not significant and it is removed from the undirected graph. 
This procedure solves the problem of false positive edges in the 
undirected graphs estimated by the glasso and by extension DCR. 

This bootstrap inferential procedure is similar to the subsam- 
pling stability selection approach of Meinshausen and Bilhlmann 
(2010). The goal is to control the familywise type I multiple test- 
ing error in a high dimensional setting by looking at the selection 
probabilities of every variable (or edge) under subsampling. In 
their framework, the data are subsampled many times and they 
choose all variables that occur in a large fraction of the resulting 
selection sets. They retain variables with a high selection proba- 
bility and remove those with low selection probabilities. The exact 



cutoff Ttthr is a tuning parameter but they notice that the results 
do not vary much for sensible choices of it. 

2.6. COMPARING GRAPHS 

Usually in an fMRI study, several subjects are scanned. In the 
multi-subject DCR (Cribben et al., 2012), we assumed the con- 
nectivity between brain regions was similar for each subject. By 
combining information across subjects, the multi-subject DCR 
averages the connectivity patterns in the group of subjects as they 
perform a set of tasks or change psychological state. However, 
in many circumstances, this may not be wise as we risk losing 
the unique patterns of activity in the individual subjects. This 
paper shows that DCR can be run on all subjects individually 
and the change points and partition specific undirected graphs 
can be compared and tested across subjects. For example, in a 
task based fMRI experiment such as an ABAB design, we could 
check whether the connectivity in the first block is the same as 
the connectivity in the third block. 

Massa et al. (2010) introduced a Likelihood Ratio test for 
testing the equality of two or more precision matrices or undi- 
rected graphs. The test originates from the test of equal covariance 
matrices for Gaussian graphical models (Anderson, 2003). It has 
the following setup. Suppose that random samples are taken from 
s multivariate normal (MVN) populations each with p variables 
Yi, . . . , Ys and the rth population has the mean vectors jjl, and 
the covariance matrix Ej. Assuming that the ijl,s vary from popu- 
lation to population, the maximum likelihood estimators of these 
are the sample mean vectors. Assume also that the mean vec- 
tors are zero and we have (/[), j = 1, . . . , f)i observations from 
Np(0, El), (y'2), j = 1, . . . , M2 observations from Np{0, E2),..., 

and iy'n), j = 1, • ■ • , «5 observations from Np(0, Es). Let Ej"^ = 
0,1, ... , Ej^^ = 0,^. Under the null hypothesis of equal precision 
matrices, Hq : Q,i = . . . = f2s, the number of unknown parame- 
ters in the common precision matrix Oq is \p{p + 1). It can be 
shown that the maximum likelihood estimator of f2o is 

s 

! = 1 

where Q,, is the ith sample precision matrbc, n, is the number of 
observations in the ith sample and n (= «i -|- . . . -|- Mj) is the total 
number of observations. The maximized log likelihood for this 
model is then 

k = K + ^n\og{\Qo\) (8) 

where K = —jnplog{2Tt) — \np and |f2ol is the determinant of 
matrix Q,q. Under the alternative hypothesis (unequal precision 
matrices), the number of unknown parameters is \sp{p + 1). The 
maximized log likelihood for this model is then 

' 1 

h=K + Y^-n,\oi{\0,\) (9) 

j= 1 
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Assuming the null hypothesis (equal precision matrices) is cor- 
rect, the LR test statistic is 



2logA = 2(/i - Zo) = «log(|S2o|) - ^«,log(|J2ol) 
'l^ol 



i= 1 



i = 1 



1^2,1 



(10) 



This test statistic has a x-squared distribution with degrees of 
freedom equal to the number of estimated parameters under the 
alternative hypothesis minus the number of estimated parameters 
under the null hypothesis. However, the test statistics's asymptotic 
X-squared distribution and its corresponding degrees of freedom 
are not valid for sparse or regularized precision matrices which is 
the emphasis of this paper. Hence, we introduce a new bootstrap 
LR test. It is set up as follows; suppose we have s subjects each 
with p variables Yi,. . . Under the null hypothesis of equal pre- 
cision matrices for all s subjects, we can bootstrap rows of data 
points across subjects. Thus, for each subject, we create a boot- 
strap data matrix Y* = Y,_ i*,. . . ,Y,; ;» by randomly selecting with 
replacement rows of data from the original data matrices Yi = 
(Yi_ 1,. . . ,Yi_ t),. . .,Ys = (Ys_ 1,. . . ,Ys, (). For each subject we esti- 
mate a precision matrix for its new bootstrap data set Y* for the 
fuU path of k values. We choose k and its corresponding preci- 
sion matrix based on the minimum BIC. We then estimate the 
test statistic in Equation (10). This procedure is repeated a large 
number of times (say, 1000) to create a bootstrap LR distribution. 

3. SIMULATIONS 

In order to assess the performance of the new DCR algorithm 
and new inference procedures, a number of simulations are 
performed. In Cribben et al. (2012) the simulation study con- 
centrated on multi-subject data but here we focus exclusively on 
single-subject data. 

The first simulation illustrates the application of the new DCR 
algorithm to identically distributed data (i.e., null data). The next 
simulations illustrate the application of the new DCR to vector 
autoregression (VAR; Zellner, 1962; Hamilton, 1995) data. For 
applications to simulated Multivariate Normal data, please refer 
to Cribben (2012). VAR is an econometric model, generalizing 
the univariate AR model, commonly used to capture the evo- 
lution and the interdependencies between multiple time series. 
More details on VAR models are included in section A.2. of the 
Appendix. VAR data are representative of the properties underly- 
ing fMRI data in that they have the property of autocorrelation 
within the individual time series (brain regions) but can have 
the possibility of non-zero cross correlation and lagged cross cor- 
relation between the time series (brain regions). For each VAR 
simulation, for regions not apart of the connected system, their 
time series are made up of i.i.d Gaussian noise indicating a lack of 
connectivity. 

The objective of each simulation is to find the times of the 
functional connectivity change points and to estimate the con- 
nectivity structure within each partition assuming that the loca- 
tion and number of change points are unknown. The stationary 



bootstrap (section 2.3) and bootstrap (section 2.5) procedures are 
used in order to perform meaningful inference procedures on the 
change points and on the edges (connections) between the brain 
regions, respectively. 

4. EXPERIMENTAL DATA 

We use the same data set here as in Cribben et al. (2012). The 
data was taken from an anxiety- inducing experiment (Lindquist 
and McKeague, 2009; Wager et al, 2009a,b). The experimen- 
tal task was a variant of a well-studied laboratory paradigm for 
eliciting social threat, in which subjects are instructed to deliver 
a speech under evaluative pressure. The design was an off-on- 
off design, with an anxiety-provoking speech preparation task 
occurring between resting periods. Participants were informed 
that they were to be given 2 min to prepare a 7 min speech, 
and that the topic would be revealed to them during scanning. 
They were informed that after the scanning session they would 
deliver the speech to a panel of expert judges, though there was 
"a small chance" they would be randomly selected not to give 
the speech. After the start of fMRI acquisition, participants rested 
for 2 min by viewing a fixation cross. At the end of this period, 
participants viewed an instruction slide for 15 s that described 
the speech topic, which was "why you are a good friend." The 
slide instructed participants to be sure to prepare enough for the 
entire 7 min period. After 2 min of silent preparation, another 
instruction screen appeared (a relief instruction, 15 s duration) 
that informed participants that they would not have to give the 
speech. An additional 2 min period of resting baseline completed 
the functional run. 

Data was acquired and preprocessed as described in previous 
work (Wager et al., 2009a). During the course of the experiment 
a series of 215 images were acquired (TR = 2 s). In order to cre- 
ate ROIs, time series were averaged across the entire region. The 
data analyzed here consists of four ROIs and heart rate for N = 23 
subjects (Figure 8A). The regions in the data set were chosen due 
to the fact that they showed a significant relationship to heart rate 
in an independent data set. The temporal resolution of the heart 
rate was 1 s compared to 2 s for the fMRI data. Hence, the heart 
rate was down-sampled by taking every other measurement. 

The new DCR approach was carried out on all 23 subjects 
individually and we were interested in testing whether stressor 
onset was associated with changes in the connectivity between 
the brain regions and heart rate. Each subject's change points 
and partition-specific undirected graphs are compared in order 
to observe whether the change points and connectivity structures 
are common across subjects. We also perform the LR test on 
the similarity of the precision matrices across subjects. With this 
information we can make inferences about whether the assump- 
tion of the same connectivity structure across subjects is valid and 
helpful. 

5. RESULTS 

For each simulation, in order to guarantee reasonable estimates 
of the sample mean vectors and sample covariance matrices, we 
set 8 = 35, that is, we set the minimum distance between change 
points to be 35 data points. This length is necessary for the single- 
subject data unlike the multi-subject data simulations in Cribben 
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et al. (2012). If 8 < 35 false positive change points occur due to 
the instability of the mean and covariance estimates. For each 
simulation we plot the BIC reduction against time in order to 
find the significant change points. In these plots, the vertical black 
lines represent the BIC reduction at the particular time point 
and the red triangles in the simulated data sets represent the 
0.975 and 0.025 empirical quantiles of the BIC reduction based 
on 1000 stationary bootstrap relications of the data for each time 
point with a non-zero BIC reduction. For the stationary boot- 
strap, we use \|/ = 1 /0.05, that is, the average block length is equal 
to 20 time points for a data set containing 100 time points. The 
new DCR algorithm stationary bootstraps the data within each 
partition after the candidate change point have been found so 
the average length of the block changes depending on the size 
of the partition. Also, for all simulations, the empirical quan- 
tile threshold for the bootstrap distribution of the edge values 
is 0.75 based on 1000 runs (stability selection; Meinshausen and 
Btihlmann, 2010). This cutoff value provided the right balance 
for including the correct edges and removing the false positive 
edges in the undirected graphs. The correct undirected graphs are 
also plotted for comparison purposes. The black and red edges 
in the undirected graphs represent positive and negative partial 
correlation, respectively. The thickness of the edges are directly 
related to the strength of the connectivity between two brain 
regions. 



5.2. SIMULATION 2 

Figure 2A depicts the splitting times plotted against BIC reduc- 
tion for a VAR simulation with four connectivity change points. 
The new DCR correctly finds all four connectivity change points. 
Figure 2B shows the true partial correlations for this VAR data 
set for the correct change points. Figures 2C,D show the undi- 
rected graphs for the partitions found in Figure 2A without and 
with inference performed on the edges using the simple boot- 
strap procedure, respectively. It is evident from the graphs that 
by performing no inference on the edges, many false positive 
edges are estimated. This simulation highlights the need for 
carrying out the bootstrap inferential procedure on the edges 
(section 2.5). 

5.3. SIMULATIONS 

This simulation is very similar to simulation 2 except in this case 
the data was generated using a VAR model with seven change 
points. Figure 3A depicts the splitting times plotted against BIC 
reduction. The new DCR correctly finds all the seven con- 
nectivity change points. Figure 3B depicts the true undirected 
graphs. Figures 3C,D show the undirected graphs for the parti- 
tions found in Figure 3A without and with inference performed 
on the edges using the bootstrap procedure, respectively. Again 
by not performing inference on the edges, false positive edges are 
estimated. 



5.1. SIMULATION 1 

As the data is i.i.d with no functional connectivity between the 
brain regions, the new DCR correctly finds no BIC reduction 
for any time point (Figure lA) and so the functional connectiv- 
ity between the brain regions remains constant throughout the 
whole time series. As the data is white noise, the new DCR cor- 
rectly finds no networks or functional connectivity between the 
brain regions (Figure IB). DCR was also applied to the same sim- 
ulated i.i.d data but with 10 spikes of magnitude 4 randomly 
added to the entire series. Again, it correctly finds no change 
points and no edges or functional connectivity between the time 
series (the figures are the same as Figures 1A,B and are not 
included). 




200 300 
Time 



FIGURE 1 I (A) Tlie splitting times plotted against BIC reduction for the i.i.d 
data set (simulation 1) using the new DCR. (B) The corresponding 
undirected graph for this data set. 



5.4. SIMULATION 4 

This simulation shows how the new DCR performs in a weaker 
signal VAR environment. The true partial correlation structure 
is shown in Figure 4C. The only difference between Partition 1 
(1-100 time points) and Partition 2 (101-200 time points) is 
that the magnitude of the partial correlation structure between 
brain regions 3, 6, and 9 is reduced while the correlation struc- 
ture between 1, 4, 8, 14 remains the same. For Partition 3 
(201-300 time points), the magnitude of the partial correlation 
structure between brain regions 1,4, 8, 14 is reduced while the 
partial correlation structure between brain regions 3, 6, and 9 
reverts back to the same correlation structure as in Partition 1. 
The connectivity structure in Partition 4 (301-400 time points) 
is i.i.d noise while Partition 5 (401-500 time points) has the 
exact same connectivity structure as Partition 1. This simula- 
tion depicts how the new DCR performs when the changes in 
connectivity are very small and in an ABA fMRI experimental 
setup. 

The new DCR correctly identifies the correct change points 
but does not deem the change point at time point 100 to be sig- 
nificant (Figure 4B). Figure 4A shows the BIC reduction plotted 
against the splits using the old DCR algorithm. Notice the dif- 
ference in magnitude of the BIC reduction and the stationary 
bootstrap confidence bounds. By recalculating the BIC reduction 
between change points a more accurate estimate is computed. The 
new DCR finds three of the four connectivity change points while 
the old DCR finds only two of the four. The new DCR method 
does a very good job of estimating the correct undirected graphs 
as well (Figure 4D). As the first change point is not deemed sig- 
nificant, the first estimated undirected graph is simply an average 
of the connectivity of the first two partitions. 
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FIGURE 2 I (A) The splitting times plotted against BIC reduction for tfie VAR 
data set with four connectivity change points (simulation 2) using the new 
DCR algorithm. The vertical black lines represent the BIC reduction at the 
particular time point and the red triangles represent the 0.975 and 0.025 



empirical quantiles of the BIC reduction based on 1000 stationary bootstrap 
relications of the data for each split. (B) The true undirected graphs for each 
partition. The undirected graphs for this data set using the new DCR and (C) 
no inference on the edges and (D) inference on the edges. 




FIGURE 3 1 (A) The splitting times plotted against BIC reduction for 
the VAR data set with seven connectivity change points (simulation 
3) using the new DCR algorithm. (B) The true undirected graphs for 



each partition. The undirected graphs for this data set using the 
new DCR and (C) no inference on the edges and (D) inference on 
the edges. 



5.5. SIMULATIONS 

This simulation is similar to Simulation 4 in that the con- 
nected brain regions are constant for the entire time course but 
the magnitude and direction of the connectivity changes over 
time. This simulation contains more change points and more 
time points. Figure 5C shows the true undirected graphs for 
each partition. The new DCR correctly identifies all the change 
points (Figure 5B) and also performs very well in identifying the 



connectivity structure for all partitions (Figure 5D). Figure 5 A 

shows the BIC reduction plotted against the splits using the old 
DCR. Notice the difference in magnitude of the BIC reduction 
and the stationary bootstrap confidence bounds. By recalculating 
the BIC reduction between change points a more accurate esti- 
mate is found. This simulation highlights the fact that the new 
DCR is very adept at locating small differences in the connectivity 
structure between two partitions for VAR data. 
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FIGURE 4 I The splitting times plotted against BIC reduction forthe VAR data set with four connectivity change points (simulation 4) using (A) the old DCR algorithm 
and (B) the new DCR algorithm. (C) The true undirected graphs for this data set and (D) the estimated undirected graphs for this data set using the new DCR. 



5.6. SIMULATIONS 

This simulation is representative of the data in the real fMRI data 
that follows in the next section. It contains five brain regions 
and 215 time points and three connectivity change points. The 
true connectivity structure for this simulation can be seen in 
Figure 6B. The new DCR correctly identifies all the change points 
(Figure 6A) and also performs very well in identifying the con- 
nectivity structure for all partitions (Figure 6C). This simulation 
demonstrates how well the new DCR algorithm performs when 
the number of brain regions is small and the number of observa- 
tions in the experimental time course is also small. 

5.7. SIMULATION 7 

Finally, in order to show the importance of the new single- 
subject DCR methodology, we present an example where the 
multi-subject DCR algorithm (Cribben et al., 2012) fails to find 
the correct connectivity change points and the correct networks 
between brain regions but the new single-subject DCR discussed 
in this work performs very well. The setup of the simulation 
is as follows — we simulate a VAR model with 500 time points 
(r = 500) and 15 ROIs (p = 15) for 10 subjects, with the first five 
subjects having positive connectivity (average correlation ~0.67) 
between ROI 8 and 15 for the first 250 time points and positive 
connectivity (average correlation ~0.7) between ROI 2 and 13 for 



the final 250 time points. On the other hand, the second five sub- 
jects having negative connectivity (average correlation ~ —0.67) 
between ROI 8 and 15 for the first 250 time points and nega- 
tive connectivity (average correlation ~ —0.7 ) between ROI 2 
and 13 for the final 250 time points. The rest of the ROIs con- 
sist of i.i.d Gaussian noise. By running the multi-subject DCR 
algorithm (Cribben et al., 2012) on this data set and aggregat- 
ing information across all 10 subjects, the two groups of subjects 
cancel each other out and no connectivity change points are 
found (Figure 7A) as well as no connectivity between ROI 8 
and 15 or ROI 2 and 13 as well (Figure 7D). In this example, 
it is not a good idea to average the information across all 10 
subjects. 

Figures 7B,C show the BIC reduction plotted against the split- 
ting time points using the new single-subject DCR algorithm for 
subject 1 (subject from group 1 which has positive connectivity 
between ROI 8 and 15 for the first 250 time points and posi- 
tive connectivity between ROI 2 and 13 for the final 250 time 
points) and for subject 6 (subject from group 2 which has negative 
connectivity between ROI 8 and 15 for the first 250 time points 
and negative connectivity between ROI 2 and 13 for the final 250 
time points), respectively. All subjects within group 1 and group 
2 behave similarly. Note that by running the new DCR on the 
individual subjects we find the correct connectivity change points 
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FIGURE 6 I (A) The splitting times plotted against BIC reduction for the VAR data set (simulation 6) using the new DCR. The corresponding (B) true undirected 
graphs and (C) estimated undirected graphs for this data set. 



and the correct undirected graphs (Figures 7E,F). This simula- 
tion illustrates how vital it is to estimate connectivity changes for 
the individual subjects and to carry out the estimation as accu- 
rately as possible. After looking at the individual subjects, we can 
then make an informed decision on whether we should com- 
bine information across subjects and run the multi-subject DCR 
algorithm. 

5.8. EXPERIMENTAL DATA 

For the fMRI data set, in order to guarantee reasonable estimates 
of the sample mean vectors and sample covariance matrices, we 
set 8 = 40, that is, we set the minimum distance between change 
points to be 40 data points. For the stationary bootstrap, we use 
\|/ = 1/0.05, that is, the average block length is equal to 20 time 
points for a data set containing 100 time points. After conclusion 



of the algorithm, we plot the BIC reduction against time in order 
to find the significant change points. Figure 8B shows the signif- 
icant change points for all subjects in the 4 ROIs and heart rate 
data set. The y-axis depicts the subject number while the x-axis 
shows the splitting times. Every subject has either two or three 
significant change points with each subject having a change point 
in the neighborhood of time point 60 ( 120 s), which corresponds 
directly to the presentation of the first visual cue specifying the 
topic of the speech and the removal of said cue at time point 67.5 
(135 s). While all subjects have a significant splitting time in the 
vicinity of the first cue, this pattern does not continue across all 
subjects for the duration of the experiment. After this, the sub- 
jects split into two groups; group 1 (subjects 1, 3, 9, 14, 15, 16, 
18, 19, 20) have only one further significant splitting time near 
time point 130 while group 2 (subjects 2, 4, 5, 6, 7, 8, 10, 11, 
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FIGURE 7 1 The splitting times plotted against BIC reduction for the 
VAR data set with one connectivity change point (simulation 7) using 
(A) the multi-subject DCR algorithm, (B) the new single-subject DCR 



(on one subject from group 1), (C) the new single-subject DCR (on 
one subject from group 2), and (D-F) their corresponding undirected 
graphs. 



12, 13, 17, 21, 22, 23) have two further significant change points 
between the first cue and the end of the experiment. At time 
point 130 (260 s), the second visual cue stating that the partici- 
pants would in fact not have to give the presentation to the expert 
panel of judges after the conclusion of the scanning session was 
revealed. 

As every subject has a significant change point in the vicinity 
of the first visual cue. Figure 9 shows the undirected graphs for 
the first partition for every subject corresponding to their individ- 
ual change points in Figure 8B. The empirical quantile threshold 
for the bootstrap distribution of the edge values is 0.75 based on 
1000 runs following evidence of stability from the simulations. 
The black and red edges in the undirected graphs represent pos- 
itive and negative partial correlation, respectively. The thickness 
of the edges are directly related to the strength of the connectivity 
between two brain regions. 

The connectivity structure in the undirected graphs for each 
subject is very different (Figure 9). The new bootstrap LR test 
discussed in section 2.6 was carried out in order to perform a 
test on the similarity of the precision matrices. The test statistic, 
A = 1860, results in a p-value of 0 rejecting the null hypothe- 
sis of equal precision matrices. However, the individual partial 
correlations shown in Figure 9 do have some overlap between 
subjects. 

For subjects in Group 1, we also compared their precision 
matrices while the participants were silently preparing their 
speeches (time points 67.5-130) and their precision matrices 
while the subjects rested after they were informed that they would 
not have to deliver the speech. In each case, we rejected the null 
hypothesis (p-value = 0) of equal precision matrices implying 
significant heterogeneity in connectivity across subjects. 

6. DISCUSSION 

DCR is an approach for splitting the experimental time course 
in a functional neuroimaging experiment into partitions based 
on functional connectivity changes between ROIs or voxels. In 
this work, we focus exclusively on the analysis of single-subject 



data, and extend the methodology in two ways: (1) we alter 
the algorithm in order to make the change point estimation 
more accurate with a small number of observations and (2) we 
perform inference on the edges in the undirected graph or con- 
nections between brain regions in order to reduce the number of 
false positive edges in the graphs. The new algorithm could be 
used for finding functional connectivity change points for multi- 
subject data as well. However, in was shown in Cribben et al. 
(2012) that the old DCR performed very well on multi-subject 
data and did not suffer from the issue of false positive edges in 
the undirected graphs to the same degree as the single-subject 
analysis. 

DCR can be applied directly to data from ROI studies or to 
temporal components obtained from a PCA or ICA analysis. It 
is very flexible in that it does not require prior knowledge of the 
nature of the experimental design and maybe particularly appro- 
priate for studies when it is not possible to replicate experimental 
manipulations within subjects. For this reason, on going collab- 
orations include applying the method to data of A-ha! or Eureka 
moments in order to observe the connectivity changes as the sub- 
ject approaches the insight moment and where the state changes 
occured. It is also being applied to resting state data, pain data 
and data from hallucination experiments. We believe the intro- 
duction of the new DCR methodology for single-subject analysis 
as well as the new bootstrap Likelihood Ratio test for sparse pre- 
cision matrices will be particularly beneficial to resting state data 
analysis as it is well known that there exists large heterogeneity 
in connectivity across subjects and between runs for the same 
subject. 

The simulation study carried out indicates that the new DCR 
method is very capable of finding small and subtle changes in 
the connectivity structure as well as robust to the inclusion of 
spikes in that data. The simulations also show that the addition 
of weaker signals do not adversely affect the results. 

In the original DCR paper (Cribben et al., 2012), the method 
was predominantly validated on multi-subject data by stacking 
subjects on top of one another in the algorithm. By combining 
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FIGURES I (A) The four-ROI and heart rate data set— the regions are: (1) VMPFC, (2) anterior mPFC, (3) Striatum/pallidum, and (4) DLPFC/IFJ. (B) The splitting 
tinnes plotted against BIC reduction for the fMRI data set (four ROIs and heart rate) for all 23 subjects using the new DCR. 
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FIGURE 9 I The undirected graphs for the first partition of all subjects in the four ROIs and heart rate data using the new DCR. The heart rate has label 5 in the 
undirected graphs. 



information across subjects, we average the connectivity pat- 
terns in the group of subjects as they perform a set of tasks 
or change psychological state. Hence, we assume that they have 
the same connectivity structure. However, group average results 
may belie the unique activity patterns of individual subjects and 
by analysing individual subjects we may be able to determine 
if an individual scan is normative with respect to a reference 
population or we may be able to understand the sources of 
intersubject variability in brain activity (Van Horn et al, 2008). 
Simulation 7 illustrates the importance of the new single-subject 
DCR algorithm. By applying the multi-subject DCR algorithm 
to this data set, the first group of subjects with positive connec- 
tivity is negated by the second group of subjects with negative 



connectivity resulting in no connectivity change points being 
found and no significant networks being found in the undirected 
graphs. The correct connectivity change points were found by 
running the new single-subject DCR algorithm on the individ- 
ual subjects. Hence, by looking at the individual subjects first 
we can make an informed decision about whether it is wise 
to run the multi-subject DCR algorithm and combine infor- 
mation across the group. This new single-subject allows us to 
perform this single-subject analysis accurately and reliably. Also, 
after performing the new DCR algorithm on the individual sub- 
jects from an fMRI study of state anxiety, we compared and 
tested the precision matrices or graphs across the individuals 
using a Likelihood Ratio test. Although the graph had some edges 
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in common, we rejected the null hypothesis of equal precision 
matrices. 

The choice of A in DCR is of particular importance as it rep- 
resents the minimum possible number of time points between 
adjacent change points. Ideally, we want to make A as large as 
possible to ensure that there is enough data to provide reliable 
estimates of the sample means and sample covariance matri- 
ces. However, the choice of A places an upper bound on the 
number of change points that can be found using our method, 
with small values allowing for more change points. Hence, we 
want to make A as small as possible to ensure that we find 
all possible change points. Ultimately, its value can be adjusted 
depending on the existence of a priori knowledge about the spac- 
ing of changes in functional connectivity. For the multi-subject 
analyses we were able to choose a A of eight time points but 
for the single-subject analyses in this paper the smallest pos- 
sible A was 30 time points. Hence, the single-subject analyses 
should only be applied to slowly changing state studies or where 
the distance between state changes was greater than 30 time 
points. 

A common question asked about DCR is how long the com- 
putation takes and how many ROIs it can handle. It is difficult 
to estimate precisely given that the alogrithm depends on sev- 
eral interconnected variables. The method's computation time 
depends on the number of time points in the data set, the num- 
ber of change points, the number of ROIs, the permutation or 
bootstrap procedures and the size of A. As the first four increase 



in number the computation time increases while as A decreases 
the computation time increases. Also, as the number of ROIs 
increase, interpretation of the undirected graphs becomes more 
difficult and the number of changes in the functional connectivity 
increases as well. The code is currently being converted to Matlab 
and paralellized in order to reduce the computational burden. We 
plan to make it available on our websites. 

The permutation and stationary bootstrap inferential pro- 
cedures for determining whether or not a change point was 
significant were reconfigured for use with the new single-subject 
DCR. In addition to inference on the change points, we also 
introduced a new bootstrap inference procedure for the edges 
or connections between brain regions in order to remedy the 
issue of excessive false positive edges in the estimation of undi- 
rected graphs after the change points had been found. The 
simulation study showed how vital this extra step is for the inter- 
pretation of the undirected graphs between each pair of change 
points. 

In conclusion, the new single-subject DCR method, intro- 
duced in this paper, will be especially helpful and informative to 
people interested in functional connectivity changes in individ- 
uals and to determine variability across subjects. We are hopeful 
this is the beginning of further work in the identification of emo- 
tion, stress or resting state transitions and aberrant connectivity 
patterns and will help our understanding of the human brain 
and neurological disorders such as depression and Alzheimer's 
disease, and the treatment of these disorders. 
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A. APPENDIX 

A.1. NEW DCR ALGORITHM 

The algorithm for performing single-subject DCR is set up as 
follows [(7-10) are the extra steps]: 

1. Consider the foil data set. Calculate the sample mean 
and covariance matrix using Equations (1, 2), respectively. 
Obtain an estimate of f2 = E^' for each value of k using 
Equation (3). 

2. Choose the value of k and hence the estimate of f2, that 
minimizes the BIC (Equation 4). 

3. Refit the model without the lasso penalty term (i.e., the li- 
constraint) keeping the zero elements in the matrix foced and 
record the minimum BIC. 

4. Partition the data into two parts: a left subset consisting of 
time points {l:y} and a right subset consisting of {y + 1:T} 
where y = A. Repeat Steps (1-3) for both data sets and sum 
the BIC scores from the two partition elements. Repeat this 
procedure for y values from A + ltoT— A+l. 

5. If the sum of the combined BIC scores for any two subsets 
is less than the BIC of the entire data set computed in (2), 
the time point with the largest decrease is selected as the 
candidate change point, thus partitioning the data into two 
segments. 

6. Apply Steps (1-5) recursively to each partition until no 
partition element can be forther split into smaller elements. 

7. Order the candidate change points once the algorithm is 
exhausted. 

8. For the first candidate change point p, partition the data into 
two parts: a left subset consisting of time points {l:p) and 
a right subset consisting of {p + l:v} where v represents the 
second candidate change point. Repeat Steps (1-3) for both 
subsets and sum the BIC scores from the two. 

9. If the sum of the combined BIC scores for both subsets is 
less than the BIC of the entire data set consisting of time 



points {l:v}, the time point remains a candidate change 
point. 

10. Repeat steps (7-9) for aU change points. 

A.2. THE VECTOR AUTOREGRESSION (VAR) MODEL 

The VAR model is one of the most used models for the analysis 
of multivariate time series models. It is used predominantly in 
the field of economics and finance to describe the dynamics 
of prices and markets. For a (« x 1) vector of time series vari- 
ables Yf = (yif, y2t, ■ ■ ■ ,ynt)^, the /"-lag vector autoregressive 
(VAR(p)) model has the form 



Y( = c -F niYt_ 1 + n2Yt_2 + . . . + npYt_p -F e,, f = l 



(Al) 



where IF,- are (nxn) coefficient matrices and €t is an (n x 1) 
unobservable mean white noise vector process with time invari- 
ant covariance matrix E (Zivot and Wang, 2003). For example, 
the equation for a univariate VAR(2) model is 
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(A2) 



yit = ci -F Jt}iyif_i -F Aiy^t-i + AiYit-i + Aiy2t-i + ^it 

yit = ci + AiYit- 1 + Aiyit-i + Aiyit-i + Aiyit-i + 



(A3) 



where covi^u, 62t) = ou for f = s; 0 otherwise. 



Frontiers in Computational Neuroscience 



www.frontiersin.org 



October 2013 | Volume 7 | Article 143 | 15 



